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SPARSE  MATRIX  TECHNIQUES  IN  TWO  MATHEMATICAL  PROGRAMIING  CODES 

by 

George  B.  Dantzlg* 

Roy  P.  Harvey t 
Robert  D.  McKnlghtt 
Stanley  S.  ^Ithf 

Introduction.  The  M3  and  M5  linear/separable  prograamlng  codes  [1]  and  [2]  were 
written  for  the  Iffif  7090/94  class  of  computer  using  the  FORTRAN  Monitor  System. 

They  solve  problems  up  to  about  400  rows.  M3  was  developed  first,  during  the 

period  1960-62,  In  part  by  RAND  Corporation  personnel  and  In  part  by  the  Standard 

Oil  of  California  (SOCAL).  It  lias  been  available  through  the  SHARE  library  for 

several  years.  MS  was  developed  later,  during  the  period  1962-64,  by  the  authors 

of  this  paper  and  Is  a  SOCAL  proprietary  code.  ! 

; 

Through  the  use  of  sparse  matrix  techniques  and  other  devices,  MS  can 
solve  a  problem  In  about  one  quarter  of  the  time  taken  by  M3.  Before  discussing 
these  Improvements,  we  list  first  some  of  the  features  that  the  codes  have  In  common: 
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Single  Precision 

Upper  Bounding  Algorithm  [3] 

Cost  Ranging  Algorithm  [4],  [5] 

Composite  Algorithm  [6] 

An  Upper-Bounded  and  a  Non-Upper-Bounded  Separable-Programming  Algorithm 
[7],  [8],  [9],  [10],  [11] 

The  ability  to  designate  free  variables  by  use  of  controls  (a  free  variable 
Is  unrestricted  as  regards  sign) 

The  ability  to  designate  frozen  variables  by  use  of  controls  (a  frozen 
variable  must  be  zero  in  the  final  solution  whether  In  or  out  of  the  basis  ) 
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Gecoff  and  Restart  Procedures 


The  ability  to  handle  multiple  objectives  and  right  hand  sides  without 
the  need  to  recalculate  the  Inverse 

Row  and  Column  Error  Calculations  as  desired 

An  Automatic  Tolerance  Regulation  Mode  of  Operation  [12] 

A  Pivot  Rejection  Algorithm 

Crashing  Procedures 

Multiple  Pricing:  The  selection  (and  possible  use)  of  several  attractive 
non-baslc  columns  during  the  prlclng-out  operation  Iti  a  simplex  Iteration 


The  most  Interesting  aspects  of  the  codes  with  respect  to  the  topic  of 
sparse  matrix  methods  are  In  connection  with  the  Inversion  techniques.  The  search 
for  more  and  more  efficient  methods  for  computer  solution  of  large  linear  programming 
problems  has  resulted  In  procedures  that  for  the  most  part  are  really  only  variants 
of  the  simplex  method.  Some  of  these  variants,  such  as  decomposition,  represent 
fairly  radical  departures;  others  differ  only  In  computational  form.  The  so-called 
revised  simplex  method  is  perhaps  an  example  of  the  latter  group  [13]. 

In  review,  recall  tliat  the  "tableau"  of  the  original  simplex  method.  Is,  in 
effect,  the  m  x  n  matrix: 

C  -  b"^a  (1) 

wliere  B  Is  the  matrix  of  basis  vectors  and  A  is  the  matrix  of  coefficients 
augmented  by  the  rlght-hand-slde. 

The  term  "revised"  of  the  "revised  simplex  method"  is  used  because  the  matrix 
multiplication,  B  ^A,  is  not  performed.  Instead,  the  Inverse  B  ^  is  maintained 
explicitly,  and  is  used  to  compute  only  those  elements  of  C  that  are  critical  to 
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Llie  simplex  algorithm.  Thus,  If  the  top  row  of  C  represents  the  relative  cost 
factors  (6j's),  we  can  compute  them  by  applying  the  top  row  of  B  ^  to  A. 

The  minimum  of  these  6j's  determines  the  pivotal  column  of  C  which  we  compute 
by  applying  B  ^  to  the  corresponding  column  In  A.  In  a  formal  sense  then,  the 


original  simplex  method  requires  us  to  recompute  an  m  x  n  matrix  each  iteration 


while  the  revised  form  requires  us  only  to  recompute  one  row  and  one  column  and  to 
update  the  right-hand-side  and  the  m  x  m  Inverse  matrix  B 


This  apparent  advantage  Is  sometimes  less  than  real,  but  the  revised 
form  offers  other  benefits  that  aid  greatly  In  the  solution  of  large  problems. 

If,  for  example,  C  Is  too  large  to  be  kept  In  the  computer's  high  speed  memory 
device,  the  revised  form  would  have  obvious  advantages.  When  high  speed  memory  Is 
no  longer  adequate  to  hold  B  the  revised  form  Is  In  trouble.  The  product  form 
of  the  Inverse  Is  a  help  In  this  case  [14]. 

As  noted,  with  the  revised  simplex  method  we  do  not  update  C;  we  only 
update  ,B  C  given  by  (1)  Is  never  computed.  With  the  revised  simplex  method, 
with  the  inverse  in  product  form,  we  do  not  update  B  ^  either.  Instead  we  compute 
a  sequence  of  elementary  matrices  rij^,  n2***»  where  an  additional  is  computed 

on  each  "re- inversion"  iteration  and  on  each  simplex  Iteration.  After  re-lnversion 
and  k  simplex  Iterations  B  ^  is  given  by  (2)  below.  It  is  never  computed. 

^  ■  Vk 

Wlien  the  price  vector  (i.e.  the  top  row  of  B  Is  to  be  computed,  the  product 
terms  of  (2)  are  premultlplled  by  p  =  (1,  0,  0,...,0)  and  the  multiplications 
executed  by  working  from  left  to  r^ght  so  that  each  multiplication  Is  that  of  a 
1  X  m  vector  by  an  elementary  m  x  m  matrix.  When  the  pivotal  column  of  C  is 
to  be  computed,  the  product  terms  of  (2)  are  post-multiplied  by  the  selected 


3 


column  of  A.  The  calculations  are  performed  from  right  to  left  so  that  each 
multiplication  Is  that  of  an  elementary  m  x  m  matrix  by  a  m  x  1  vector. 

Host  problems  require  so  many  Iterations  that  the  product  form  becomes 
Impractical  unless  the  product  Is  periodically  reduced  to  Its  fewest  factors. 
Unless  there  happen  to  be  unit  basic  columns,  this  minimum  number  Is  m  as  shown 
by  (2).  The  usual  practice  Is  to  "throw  away"  all  the  n's  and  to  reconstruct 
m  new  ones  from  the  current  basis  B.  In  LP  parlance  this  process  Is  called 
"re-lnveralon"  or  simply,  "Inversion". 

If  the  usual  simplex  rules  are  adhered  to  and  we  neglect  the  occurrence 
of  "ties",  the  n's  outside  the  parentheses  In  (2)  are  determined  uniquely. 

These  are  the  n'a  produced  during  the  simplex  Iterations  and  their  values 
result  from  the  simplex  pivot  criteria.  The  n's  Inside  the  parentheses  In  (2) 
are  the  ones  produced  during  Inversion  where  different  criteria  may  be  applied. 

The  most  successful  criteria  seem  to  be  those  that  seek  two  major  goals;  In  their 
purest  form  they  are: 


(a)  Numerical  accuracy 

(b)  Few  non-zero  elements 

We  will  return  to  (a)  later.  In  order  to  appreciate  (b).  It  must  be  noted  that 
most  large  LP  matrices  have  a  low  percentage  of  non-zero  elements.  (A  density 
of  less  than  5X  Is  not  unusual.)  If  the  n's  can  be  kept  relatively  sparse  they 
will  not  only  Imply  less  arithmetic,  but  It  will  also  be  possible  to  hold  them 
In  less  computer  store  by  the  use  of  a  compact  format  containing  only  the  non¬ 
zero's  (In  the  relevant  column)  and  their  row  Indices. 
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As  far  as  the  authors  know  there  Is  no  known  theoretical  way  (except 


perhaps  exhaustive  search  of  pivots)  for  finding  the  product  form  of  the  Inverse 
that  has  the  minimum  number  of  non-zero  elements.  Some  very  simple  procedures, 
however,  In  practice  have  shown  quite  remarkable  Improvements  over  those  where 
no  attention  to  (b)  Is  given.  Selecting  the  pivot  columns  of  B  In  the  order 
of  Increasing  density  often  has  a  profound  effect  on  the  percentage  of  non-zeros 
In  the  Inverse  product  form.  (In  practice,  the  entire  matrix  A  may  be  pre¬ 
ordered  by  density  which  makes  the  ordering  of  any  basis  B  automatic.) 

I 

One  of  the  measures  we  have  used  In  connection  with  criteria  (b)  shown 
above  is.  the  percent  Increase  In  non-zero  elements  In  the  representation  of  the 
inverse  compared  with  the  number  of  non-zero  elements  in  the  basis  matrix  before 
Inversion. 

M3  uses  the  product  fcrm  of  the  Inverse  and  early  versions  selected  pivot 
columns  during  Inversion  by  considering  basis  columns  In  the  order  they  were 
originally  given  and  the  pivot  row  was  selected  by  taking  the  largest  pivot  In 
absolute  value  on  the  remaining  rows.  On  typical  problems  In  the  2-5%  density 
class.  It  was  found  that  often  there  was  an  Increase  of  about  5-10  fold  In  non¬ 
zero  elements  after  Inversion.  Using  Just  the  simple  technique  of  presenting  the 
columns  by  density  (actually  only  a  partial  sort  on  large  problems)  It  was  found 
that  this  Increase  In  non-zeros  could  be  cut  by  a  factor  of  about  two. 

A  related  Idea  Is  the  so-called  "merit"  sort  which  orders  the  columns 
by  Increasing  merits.  The  total  merit  In  a  column  Is  the  sum  of  merits  assigned 
to  each  non-zero  element  In  the  column.  The  merit  of  a  non-zero  element  In  row  1 
Is  the  number  of  non-zero  elements  In  row  1.  Some  other  techniques  of  this  kind 
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are  described  In  Larsen  [IS]  and  Tewarson  [16].  The  compact  Inversion  technique 
used  In  MS  was  Influenced  by  the  above  experience  with  various  criteria.  MS  uses 
the  "Elimination  Form  of  the  Inverse",  see  Markowitz  [17].  The  notable  thing 
about  this  evolution  of  computational  Improvements  Is  that  each  one  Introduces 
something  additional  to  be  done  Implicitly  rather  than  expllcitl} :  the  revised 
method  does  not  explicitly  compute  C,  the  product  form  method  dees  not  explicitly 
compute  B  the  Markowitz'  method  does  not  explicitly  compute  the  elements 
of  In  rows  previously  used  for  pivoting.  ' 

It  has  long  been  known  that  triangularlzatlon  followed  by  back  solution 
requires  less  calculation  than  full  dlagonalization  In  lOOX  dense  systems.  Our 
experience  indicates  that  the  same  holds  for  sparse  systems.  Application  of  the 
triangularlzatlon  approach  to  a  staircase  structured  basis  Is  found  In  [18].  In 
[19]  Bartels  discusses  numerical  accuracy  of  elimination  procedures. 


Elements  on  rows  already  used  for  pivoting  are  not  eliminated  In  M5  but 
are  stored  as  a  matrix  T  which  Is  triangular  with  respect  to  the  pivot  diagonal. 

The  transformation  matrices  are  similar  in  structure  to  a  product  form 

transformation  matrix  with  the  exception  that  there  are  no  elements  above  the 
diagonal,  the  diagonal  being  defined  from  left  to  right  from  successive  pivot 
positions.  Ignoring  a  permutation  matrix  which  in  practice  is  effected  by  maintaining 
the  indices  of  pivot  positions,  the  following  relationships  hold; 


E  E  B 

m  m-1  1 


b“^  -  t“^e  E  ,...E.E, 

m  m-1  2  1 


As  with  the  product  form,  after  k  simplex  iterations,  we  will  have: 


®  “  \+k\+k-l'”^nri-l^’^  ^m®m-l’**®l^ 
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The  device  of  course  Is  not  to  compute  T  ^  but  only  T.  In  ’M5  the 
elements  of  the  triangular  matrix  T  are  stored  by  columns.  Vfhen  the  top  row 
of  Is  desired,  we  first  compute  from  left  to  right: 


^1  ■  ^Vk\+k-l***Vl 

then 


Is  found  by  "backsolvlng"  the  triangular  system  P2T  ■  and  finally,  we  compute 
the  product,  again  from  left  to  right: 


PnE  E  . • a • E, 
L  m  m-l  1 


The  selected  column  of  C  Is  calculated  in  an  analogous  right  to  left 
manner.  The  matrices  T  and  E^  are  not  unique  so  that  there  exists  the 
opportunity  to  seek  goals  (a)  and  (b).  With  this  In  mind,  we  considered  several 
pivot  selection  rules. 


Pivot  Selection 
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In  our  Initial  investigations  five  pivot  selection  algorithms  were 
considered,  which  are  described  below.  A  comparison  between  them  was  made  by 
performing  symbolic  trlangularizatlon  on  several  bases  taken  from  actual  applied 
problems  and  keeping  track  of  the  position  of  non-zero  elements.  This  symbolic 
trlangularizatlon  was  carried  out  on  the  computer.  The  proliferation  of  non-zero 
elements  was  all  that  was  recorded.  The  arithmetic  of  the  reduction  procedure 
was  not  performed.  One  consequence  of  this  Is  that  any  cancellation  on  non-zero 
elements  giving  rise  to  a  zero  In  any  particular  spot  was  assumed  not  to  occur. 
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The  five  pivot-selection  algorithixs  considered  have  one  characteristic 
In  common.  They  are  all  dynamic  In  the  sense  that  at  each  stage,  the  next 
pivot  is  chosen  In  the  light  of  the  updated  matrix.  In  this  way,  the  criteria 
differ  from  some  of  the  techniques  mentioned  earlier.  At  each  plvot-selectlon 
step,  the  choice  depends  only  on  properties  of  that  part  of  the  updated  matrix 
made  up  of  rows  and  columns  not  yet  made  use  of  for  pivoting  and,  of  course, 
the  pivot  choice  Is  made  from  this  submatrix.  The  number  of  non-zero  elements 
In  each  such  partial  row  and  column  Is  recorded  and  kept  currient  at  each  pivot 
step. 


The  flvbt-selection  criteria  considered  are  now  described.  Of  all  the 
positions  In  the  as  yet  unplvoted  portion  of  the  updated  matrix  which  contain 
a  non-zero  element: 


1)  Select  as  next  pivot  row,  the  row  with  the 
minimum  number  of  non-zero  entries.  In  the 
case  of  ties,  take  the  first  such  row.  As 
the  pivot  column,  take  the  first  column 
which  contains  a  non-zero  entry  In  this 
pivot  row. 

2)  Select  next  pivot  row  as  In  1).  Of  all  columns 
which  have  a  non-zero  entry  In  the  pivot  row, 
select  the  first  column  which  contains  the 
minimum  number  of  non-zero  entries. 

3)  Consider  the  subset  of  rows  {R}  which  have 
the  minimum  number  of  non-zero  entries.  There 
will  be  only  one  such  row,  unless  ties  occur. 

Of  all  the  columns  which  have  at  least  one  non¬ 
zero  entry  on  a  row  of  {R},  select  the  first 
column  which  contains  the  minimum  number  of  non¬ 
zero  entries.  This  is  the  pivot  column.  The 
first  row  of  {R}  containing  a  non-zero  entry 
In  the  pivot  column  Is  taken  as  the  pivot  row. 

4)  Select  that  element  as  pivot,  such  that  the 
product  of  the  number  of  other  non-zero  elements 
In  Its  row,  times  the  number  of  other  non-zeros 
In  Its  column.  Is  a  minimum.  This  criterion  is 
due  to  H.  Markowitz  [17]. 
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5)  Take  that  position  containing  a  non-zero  entry 
as  pivot  which  creates  the  fewest ' additional 
non-zero  entries,  when  carrying  out  the 
reduction  step. 

6)  to  ID).  Apply  the  same  rules  as  above  to  the 
transposed  matrix.  Note  that  9)  could  lead  to 

a  different  choice  of  pivot  than  4  because  of  the 
tie-breaking  rules;  similarly  10) differs  from  5) 


In  each  case,  the  total  number  of  non-zero  entries  in  the  and  T 

I 

matrices  is  compared  with  the  number  of  non-zeros  in  the  original  matrix  B  of 
basis  vectors. 


The  results  for  two  matrices  are  shown  in  Table  1.  They  are  typical 
of  results  that  we  have  observed  on  other  matrices  which  are  not  presented  here 
because  of  incomplete  information  regarding  the  source  of  the  matrices  or 
incomplete  analysis  on  our  part.  In  our  opinion,  it  would  be  worthwhile  for  some 
one  to  collect  a  wide  variety  of  matrices  used  in  LP  applications  (not  randomly 
generated  ones)  to  verify  whether  or  not  these  results  are  typical. 
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Matrix  I 


Matrix  n 


SIZE 

62  X  62 

216  X  216 

\ 

NUMBER  OF  NON-ZERO 
ENTRIES 

348 

‘ 

1400 

DENSITY 

9% 

3X 

ALGORITHMS 

NON-ZERO  ENTRIES 
IN  E^  AND  T 

X 

INCREASE 

NON¬ 

IN 

t 

•ZERO  ENTRIES 
AND  T 

X 

INCREASE 

1)  Minimum  Row  - 
Column  with  1st  non 
zero  entry. 

378 

8.6 

1456 

4.0 

2)  Minimum  Row 
Minimum  Column 

364 

4.6 

1485 

6.1 

3)  Minimum  Rows 
Minimum  Column 

366 

5.2 

' 

' 

1424 

1.7 

4)  Minimum  Product 
[Markowitz] 

378 

8.6 

1492 

6.6 

5)  Minimum  Increase 

420 

20.7 

1553 

10.9 

6)  Same  as  1  on 
transpose 

408 

17.2 

1654 

18.1 

7)  Same  as  2  on 
transpose 

361 

3.7 

1493 

6.6 

8)  Same  as  3  on 
transpose 

361 

3.7 

1461 

4.4 

9).  Same  as  4  on 
transpose 

381 

9.5 

1514 

8.1 

10)  Same  as  5  on 
transpose 

366 

5.2 

1523 

8.8 

TABLE  1.  PIVOT  SELECTION  CRITERIA 
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Our  Intuition  would  have  probably  led  us  to  algorithm  5),  which 
did  not  show  up  very  favorably  In  our  admittedly  few  experiments.  We  chose  to 
adopt  algorithm  8)  for  the  LSP  code  MS. 

It  Is  Interesting  to  note  that  our  experiments  indicated  the  above 
technique  performs  much  less  favorably  on  sparse  matrices  whose  non-zeros  have 
been  generated  In  a  random  manner.  Symbolic  matrices  of  5Z  density  and  of  orders 

I 

50,  200  and  400  were  generated  using  the  pseudo  random  number  generator  based  on 

X. ,,  ■  aX,  (mod  2^^)  where  a  ■  X  ■  5^^ 

1+1  1  o 

The  location  of  a  non-zero  element  In  the  matrix  was  determined  by 
randomly  selecting  Its  row  and  Its  column.  To  Insure  non-singularity,  symbolic 
non-zero  elements  were  first  located  along  the  main  diagonal. 

The  results  using  algorithm  8)  are  as  follows: 


Random  Random  Random 
Matrix  1  Matrix  2  Matrix  3  Matrix  I  Matrix  H 


SIZE 

50 

200 

400 

62 

216 

DENSITY 

5% 

5Z 

5% 

9% 

3% 

INCREASE  OF  DENSITY  AFTER 
TRIANGULARIZATION  USING 
ALGORITHM  8) 

13% 

333* 

480% 

3.7% 

4.4% 

These  results  are  evidence  that  general  conclusions  should  not  be  drawn  from  linear 
programming  algorithm  experiments  which  have  been  performed  using  matrices  whose 
elements  have  been  generated  in  a  random  manner. 
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The  exact  foztn  uf  the  algorithm  Implemented  Is  given  below.  For 
numerical  accuracy,  the  above  procedure  Is  modified  If  the  pivot  chosen  Is 
less  In  absolute  value  than  a  specified  tolerance  value  (see  step  3).  The 
test  and  procedure  In  the  case  of  an  Ill-conditioned  basis  Is  omitted. 


STEPS  IN  INVERSION  PROCEDURE 


1  Count  and  store  number  of  non-zeros  In  basis  matrix  by 
column  and  row. 


2  Restrict  the  candidates  for  pivot  to  non-zero  elements 
In  columns  and  rows  not  previously  used  for  pivot  whose 
non-zero  column-counts  are  minimal. 

Further  restrict  the  candidates  to  those  whose  absolute 
values  are  ^  a  given  tolerance  value;  If  none  then 
those  whose  absolute  values  are  maximal. 

Further  restrict  the  candidates  to  those  being  In  rows 
whose  non-zero  row-counts  are  minimal. 

Among  the  latter  select  as  pivot  the  one  with  the  lowest 
column-index  and  then  the  lowest  row-lnaex. 


3  Generate  transformation  column  E^  and  corresponding 
column  of  T. 

4  Update  residual  matrix,  adjusting  the  non-zero  row  and 
column  counts  of  elements  In  non-plvoted  rows  and  columns. 

5  Repeat  steps  2  through  5  until  all  columns  have  been  used 
for  pivoting. 


The  Implementation  ofthls  algorithm  also  presented  us  with  some  Interesting 
systems  problems  In  data  handling.  We  choce  to  restrict  the  application 
of  the  code  to  problems  where  Inversions  could  be  carried  out  entirely  in  the  high 
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speed  storage  of  the  computer.  As  noted  earlier  both  M3  and  MS  are  used  for 
problems  up  to  about  400  rows.  With  the  density  and  size  problems  encountered 
MS  has  carried  the  elimination  form  of  Inverse  successfully  in  core. 

Some  comparison  times  of  M3  vs  MS  are  given  in  Table  2. 


PROBLEM  SIZE 


MS  as  Z 
of  M3 


99 

16S 

Complete  Problem 

67  secs 

42  secs 

63Z 

166 

744 

Average  Inversion 

21  secs 

10  secs 

48Z 

Average  Iteration 

3.1  secs 

1.0  secs 

32Z 

264 

1118 

Average  Inversion 

100  secs 

2S  secs 

2SZ 

Average  Iteration 

4.7  secs 

1.2  secs 

26Z 
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